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WHOLE NUMBER STRAPDOWN COMPUTATIONS 


ABSTRACT 

An inertial navigation system employing a gimballess inertial measurement 
unit requires an analytical transformation of the vehicle co-ordinate system 
into the inertial co-ordinate system. An algorithm is developed for maintaining 
an up-to-date transformation matrix in a general purpose whole number computer. 

A method of implementing the algorithm in the Apollo Guidance Computer (AGC) 

is described. The performance of the algorithm, the effects of flight 

profile parameters upon the accuracy of the algorithm, and the effects of 

certain equipment constraints are detailed in the results of computer simulations. 

Extensive computer simulations were conducted to verify the validity of the 

algorithm; while conclusions about navigation computer design were drawn 

from the simulation results, raw simulation data is included for individual 

interpretation. For comparative purposes, the results of simulation of 

a digital differential analyzer (DDA) are included. It is concluded that 

for at least certain missions, general purpose computers can be built 

which will perform the strapdown computation with sufficient accuracy 

and which will not significantly detract from the other tasks required 

of the general purpose computer by doing these tasks fast enough. 


by J. C. Pennypacker 
February, 1966 


3 




TABLE OF CONTENTS 


Section Page 

I Introduction 7 

II The Cosine Matrix 9 

III The Basic Algorithm 11 

3.1 Derivation 11 

3.2 Interpretation of Algorithm 20 

3.3 Error Accumulation 24 

3.4 Timing Considerations 25 

IV Digital Differential Analyzer (DDA) 27 

V Computer Simulation 31 

5.1 Goals of Simulation 31 

5.2 Determination of True C Matrix 32 

5.3 Characteristics of Simulated Flight Profiles 33 

5.4 Effect of Vehicle Rotation 35 

5.5 The Basic Profile 37 

5.6 Non- Synchronous Accelerations 43 

5.7 Extensions of the Basic Algorithm 46 

5.7.1 Reduction of Sampling Time Interval 46 

5.7.2 Fourth Order N Matrix 49 

5.7.3 Interrupted Sampling Time Interval 50 

5.8 Practical Considerations 56 

5.8.1 Computer Word Length 56 

5.8.2 Gyro Limitations 57 

5.9 LEM Profile 58 

5.10 Comparative Data for the Basic Profile 67 

VI Matrix Orthogonality 85 

VII Conclusions 89 


5 



I. INTRODUCTION 


There is at the present time considerable interest among designers 
of inertial navigation systems in the utilization of inertial measurement 
units which are mounted directly to the vehicle frame; the resulting con- 
figuration is a gimballess inertial measurement unit (GIMU) , commonly referred 
to as a "strapdown" system. Such a configuration requires that analytic 
methods rather than the conventional physical gimbals be employed to isolate 
the vehicle co-ordinate axes from the inertial co-ordinate system in which 
navigation and guidance of the vehicle are performed. There are at least 
two basic methods of implementing the required analytic functions: the 
more generally accepted approach is to use a digital differential analyzer 
(DDA), the other approach is to use a general purpose whole number computer. 

The desirability of the latter method becomes pronounced in those systems 
for which a general purpose computer is required to perform functions other 
than those required for navigation; in such a system, the hardware configuration 
need not include an extra processor - specifically the DDA - for navigation. 

The primary question in using a general purpose computer centers around 
the algorithms used for updating the transformation matrix. For the general 
purpose computer approach to be practical, the computer must spend only a 
small fraction (less than 10%) of its time in the strapdown task otherwise 
performed by the DDA. The time spent by the general purpose computer is 
a function of both computer speed and the updating algorithm utilized. 

The DDA algorithms are ill suited for implementation in a general purpose 
computer and the question thus arises as to whether the whole number algorithms 
of the class proposed by A. Hopkins will give sufficiently precise results 
without requiring excessive computation times. While the advantages, dis- 
advantages and capabilities of the DDA are generally understood, such insight 
into the performance of a general purpose computer operating in conjunction 
with a strapdown navigation system is lacking. 


^ Albert Hopkins, Digital Development Report #5, Updating a Cosine Matrix 

in a Whole Number Computer , MIT Instrumentation Laboratory, August 12, 1964. 
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This paper presents one algorithm which could be used to perform 
navigation functions on a whole number general purpose digital computer; 
the results of extensive computer simulation of this algorithm are also 
included. Because of the current interests of the author, the study under- 
taken is oriented towards the Apollo mission; of specific interest is the 
feasibility of utilizing the Apollo Guidance Computer (AGC) to perform the 
navigation functions of the Lunar Excursion Module (LEM). The scope of 
this study is restricted to one portion only of the general navigation 
problem, that of maintaining an accurate and timely direction cosine matrix. 
The vehicle containing the strapdown system is assumed to be a spacecraft 
of the LEM type; this assumption is fundamental to the characteristics 
of the algorithm and simulations presented herein. 
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II. THE COSINE MATRIX 


In order to perforin the navigation and guidance computations in a 
fixed co-ordinate system, it is necessary first to resolve the accelerations 
measured by the accelerometers in the spacecraft (body) co-ordinate system 
into components in the fixed co-ordinate system. For a fixed co-ordinate 
system F and a body co-ordinate system B, the transformation of acceleration 
from the body system to the fixed system is given by the following equation: 


\ = [c] ^ 


( 1 ) 


Ap = Acceleration vector resolved into the 


fixed co-ordinate system. 


[c] = Transformation matrix. 

Ag = Acceleration vector resolved into the body co-ordinate system. 


The transformation matrix is a matrix composed of the direction cosines 
of the angles between the axes of the two co-ordinate systems; thus, the 
elements of [c] are given by the following: 


C. . 
ij 




( 2 ) 


u 


( 


) 


a unit vector in the direction of the co-ordinate system 
indicated by the subscript. 


The matrix [c], which is dependent only upon the attitude of the vehicle, 
must be precisely known at the time accelerations occur in order to determine 
the position in inertial space of the spacecraft. The analytical determination 
of the C matrix is the basic difficulty encountered in the strapdown configuration. 
As the vehicle rotates, the matrix [C] changes; thus in general the velocity 
of the spacecraft in the fixed co-ordinate system is given by: 

t 

V F (t) = J [C(t)] Ag (t) dt 


0 


( 3 ) 



The inertial position of the spacecraft is determined from a further integration 

of equation (3). In order to determine an expression for the change of 

[C(t)] as the vehicle rotates, let the vehicle rotate with respect to the 

fixed co-ordinate system with an angular velocity fi (t). Then from equation (2) 

FB 


- c. .(t) = C. .(t) 

dt 1J 


Fi 


Bj 


— > 
+ u 


Fi 


Bj 


(4a) 


= U Fi • (n FB (t) X U Bj ) + ° 

Evaluating the vector equation and writing in index form yields: 


C ij (t) = U Fi ‘ [ -W t)u Bi + “ Q FBi (t)u Bk ] 


(4b) 


From equation (2) this can be rewritten as: 




th 


Letting x, y and z represent the i, j and k axis of the spacecraft 
respectively, equation (4c) can be expressed as: 


[C(t)] = [C(t)][fi(t)] 


( 5 ) 


where 


[fi(t)] = 


0 -^(t) ® y (t) 

0) z (t) 0 -^(t) 

-® y (t) ® x (t) 0 
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III. THE BASIC ALGORITHM 


3.1 Derivation 

( 2 ) 

A. Hopkins has described a method of approximating the solution 
to equation (5) utilizing a general purpose digital computer. Because this 
approximation provides the basis for the computer simulation, the remainder 
of this section presents the algorithm originally described by Hopkins. 

Define a matrix [M(T) ] which is a function of the cu' s , of their derivatives, 
and of a sampling time interval T. The matrix [M(T) ] relates the C matrix 
at the end of the sampling time interval T to its value at the beginning 
of the time interval. This relationship is defined thus: 

[C (T) ] = [C (0) ] [M(T) ] (6) 


Knowledge of M(T) enables one to calculate the current value of [C(t)] 
by a recursive process. Owing to limitations of GIMU instruments, however, 
[M(T)] can only be approximated. Previous approaches have emphasized the 
use of digital differential analyzers (DDA's) in order to achieve maximum 
precision with a small computer. Utilization of a general purpose digital 
computer such as the Apollo Guidance Computer (AGC) requires a substantially 
different approach: a large time interval T with a sophisticated approximation 
to [M(T)] instead of the DDA's short interval and skeletal approximation 
to [M(T)]. The fundamental question associated with large intervals T centers 
around the uncertainties as to the order in which rotations occur, and the 
inaccuracies which result from these uncertainties. 

The data from which [M(T) ] can be approximated by a spacecraft navigation 
computer is a quantized representation of angle changes as detected by the 
body-mounted gyros. It is here assumed that these angle changes are known 
precisely; the effect of introducing imperfect gyros into the system 
is described in a later section. To express [M(T) ] in terms of the spacecraft 
angle changes (denoted 0^, 0^, 0^) , [C (t) ] is expressed as a function of 
[C (0)3. The Taylor series expansion of element C„ is: 


( 2 ) 


Ibid. 
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( 7 ) 


• JL r\ *0 n ••• 

C..(t) = C. .(0) + t C. .(0) + - t C..(0) + - t C..(0) + . . . 

1J 1J 1J 2 1J 6 1J 

At this time it is convenient to rewrite equation (5) in index form: 


c il(t) = C. 2 (t)05 z (t) - C t3 (t)cn y (t) 

C i2 ( t ) = - C i]L (t>® (t) + C. 3 (t)co x (t) (8) 

C i 3 (t) = c il (t)03 1 (t) - c i2 (t)<n x (t) 


The expressions of equation (8) cai be used to replace the C„ (t) term of 
equation (7) with undifferentiated terms. Equation 8 can also be differentiated 
to give C „ ' s in terms of C„'s. Further substitutions utilizing equation (8) 

yield expressions for the Ch , ' s in terms of the C..'s alone. These expressions 

# ® J 3 

may be substituted for the C term of equation (7). For example: 


c il(t) = c i2 ( t >“ z ( t ) + a) z < t ) C i 2 (t) ' C i3 (t)<D y (t) “ 

= c i2 <t)“ z ( t ) - °i (t) c il (t) + C0 x (t)cD z (t)C i3 (t) 

- a5 y (t)C i3 (t) - o? (t)C il (t) + ^ x (t)03 y (t)C i2 (t) 

= E-“y (t) - <4 (t)] C il (t) + [® x (t)a, (t) + co z (t)] C. 2 (t) 

+ [03 x (t)03 z (t) - a> y (t)] C. 3 (t) (9) 


Continuing in this manner, one can obtain expressions for the time derivatives 

of each C. . in terms of all the C. ,'s. Since these expressions contain cd' s 
tj ij 

and their derivatives, they will be of the general form: 


d > V 0 ■ + f l j k2 [ ‘“< t > ]C i2 (t > + £ 1Jk3 t»(t)] c 13 (t) 


(10a) 
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Specifically, at time t = 0, equation (10a) becomes: 


d 

dt 


— V°> ' £ ijkl [a) <°> ]C il<°> + ]C i2 C°> + f l;|k3 [ “< 0 > ]C 13 (0 > 


When equation (10b) is substituted for the differentiated terms of 
equation (7), one obtains: 


C (t) = Z f [co(0) ] - C (0) 

1J k=0 1Jlci kl xL 


+ z f. jk2 MO)] c. 2 (0) 


+ S f MO)] - c (0) 

k=0 J J k'. J 


( 11 ) 


Comparison with equation (6) shows that at time t = T, the elements of 
[M(T)] are given by the infinite series in equation (11). Elimination of 
the redundant subscripts in equation (11) leads to the expression: 


M (T) = Z f M0)] - 
J k=0 Jk k! 


( 12 ) 


which is recognized as the Taylor series expansion of M. . (T) where f is 


th 


ij 


ijk 


the k derivative of M. .(0). 

It has been shown that the elements of [M(T)] can be approximated by 
a Taylor series whose terms are obtained from differentiation of equation (8); 
a list of these terms is given in Table 1. There remains to be shown how 
these terms can be expressed in terms of the spacecraft angle changes during 
the time interval T. 

til 

Letting the change of the spacecraft angle about the i axis be denoted 
by 9^, the first step is to use the Taylor series expansion to relate the 
9's to the respective co' s . According to the definition of 0^(T): 


(10b) 
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T 

e.(T) = I co i (t)dt 
o 

T 2 

/ • L • • 

[no. (0) + tax (0) + - ax (0) + . . . ]dt 

!L 1 1 

0 



= Tcd. (0) + - cd. (0) + - cd. (0) + o . . 

3 - O 1 H 


(13) 


It is evident that terms of equation (13) appear also in the Taylor expansion 
of some of the M. ,'s. For example, from equation (13) and from Table 1: 
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Order 



15 


All CD's are evaluated at time 0, the beginning of the sampling interval T. 


Order M 12 (T) M^tf) M 32 (T) 
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e (T) = To) (0) + - CD (0) + - CD (0) + . . 


M 13 (T) = Tco y (0) + - [cd x (0)cd z ( 0) + CD y (0)] 


+ - [03 (0) + 03 (0)03 (0) + 2cd (0)03 (0) - 03 (0)o3 (0)] + . . 

ft y x z x z y 

Comparison shows that 9 (T) approximates M- ^ (T) with an error function whose 
leading terms (for T < 1) are: 


03 (0)o3 (0) + - [03 (0)O3 (0) + 203 (0)o3 (0) - 03 “ (0)o3 (0)] (14) 

xz /- x z xz y 


An improved approximation to M^(T) is obtained by expressing the first error 
term of expression (14) using the product ® x O)® z (T)» From equation (13): 


9 (T) (T) = T cd (0)cd (0) + - [cd i (0 )cd (0) + CD (0)co (0)] + . . . (15) 

xz; xz 2 x ^ xz 

Utilizing Table 1, equation (13) and expression (14), one can approximate 
e x (T)0 z (T) 

M.„(T) by 9 (T) + — — with an error function whose leading terms are 

y 2 


[CD (0)CD (0) - CD , (0)05 (0) - 2 CD (0)CD . (0) ] 

Z X X Z y 


Table 2 gives a number of functions of 0. which are used in the 

1 

formulation of yet better approximations to the (T) . 
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2 3 4 

T • T .. T ... 

0 . = TO). + - CD. 4- - CD. + - CD. + . . 

1 1 2 1 6 1 24 1 ‘ 


2 2 2 3 4 ^'2 ^ " 

0 = T CD. + T CD. CD. + T (-CD. + - CD. CD.) + . . . 

1 1 11 . 1 „ 1 1 

4 3 


2 T . . 4 1 * ’ 

0 0 = T CD. CD. + - (CD.CD. + CD. CD. ) + T (- CD. CD. + 

lj i j 2 l J j i ^ l j 


- T 4 

1 2 . 1 3.. T 

0 = / CD. (t)dt = -Tcd. + - T CD. - - T J CD. + - 

-1 Q / 1 1 2 1 6 1 24 


T 1 . 

2 1 * ' 4 • 

0.0 _ = -x CD.CD. + - (CD.CD. - CD.CD.) + T (- CD.CD. 


1 -J 1 J 


13 J 1 


1 J 


3 

0 0 - 0 0 . = T (CD.CD. - CD.CD. ) + . . 

1-3 3-1 31 31 


2 3 2 4 2 ' 2 * 

00 = T CD CD. + T (- CD CD. + CD.CD. + CD.CD. .CD. 

i 1 2 1 1 1 1 1+ i i + i 


1 .0 1 

-CD.CD. + - CD.CD. ) + . . . 

6 1 J 6 J L 



1 1 .. 

- - CD.CD. - - CD.CD.) + . . „ 

6 1 J 6 1 J 


+ CD.CD. CD. ) 

1 1+2 i+2 


NOTE: All cd's are evaluated at time 0, the beginning of the sampling interval T. 

?222 
CD = CD + CD + CD 

x y z 

0 . is the negative of the angle change in the preceding interval. 


TABLE 2. AUXILIARY FUNCTIONS OF 0 
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It is now convenient to define a matrix [N(9)] which approximates 
[M(T)] with an error [E(T)]; that is: 

[M(T) ] = [N (9) ] + [E (T) ] (16) 

The number of possible forms of [N(9)] is of course unlimited and no procedure 
is given here for deriving optimized approximations. Table 3 shows three 
N matrices which are equivalent to Taylor expansions of [M(T)] to the first, 
second and third order terms respectively. The N matrices are written in 
terms of body angles; the leading terms of the corresponding error matrices 
[E(T)] are expressed as functions of the cd' s and their derivatives. The 
process of updating the N matrices of Table 3 at regular sampling time 
intervals T constitutes the basic algorithm; modifications to this basic 
algorithm are discussed later. 

3.2 Interpretation of Algorithm 

The algorithm presented in the preceding section was developed from 
a purely mathematical basis with no physical interpretation of the algorithm 
included. The transformation matrix can be visualized as a vector originating 
at the center of the unit sphere and terminating on the surface of the unit 
sphere. Rotation of the vehicle employing the strapdown system corresponds 
to tracing a path on the surface of the unit sphere. The N matrix vector, 
which approximates the true transformation vector, is updated only at discrete 
time intervals. Because the vector addition of small angle changes is 
an ordered process, the N matrix vector which is updated based upon the algebraic 
sum of angle changes during a sampling time interval is accurate only within 
some cone of error. To reduce the size of this cone of error, the position 
of the N matrix vector is extrapolated not only on the basis of historical 
velocity, but also on the basis of changes in velocity. This extrapolation 
is evidenced by the inclusion in the third order update formula of angle 
changes over two successive sampling time intervals. The physical assumption 
which is thus being made in the basic algorithm is that changes in angular 
positions during successive sampling intervals caused by acceleration are 
small compared to changes in angular position caused by current rotation. 
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NOTE: All a).' s are evaluated at time 0, the beginning of the sampling interval T. 
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TABLE 3 (cont 



3.3 Error Accumulation 

The accumulation of attitude error is a complicated process, and no 
manageable analytic description has been developed. However, a crude upper 
bound of the accumulated error can be calculated by using the assumption 
that the absolute value of the error is the sum of the absolute values 
of the errors at each update calculation. As an example of this calculation. 


consider the error terms of the matrix [N^CT)]. Since: 

[C (T) ] = [C (0) ] [M(T) ] (17) 

and 

[M(T>] = [N 3 (T)] + [E 3 (T)] (18) 

it follows that: 

[C (T) ] = [C(0)][N 3 (T)] + [C(0)][E 3 (T)] (19) 

where the second term is the error resulting from the approximation 
formula N 3> Let this error be denoted by [D(T)], i.e., 

[D(T)] = [C(0)][E 3 (T)] (20) 

Referring to Table 3 we can write: 

D il ( T) = C n e n (T) + C. 2 e 21 (T) + C. 3 e 31 (T) + 0( T 5 ) (21) 


where the e's are the elements of [E 3 (T)]. An upper bound to equation (21) 
can be obtained by substituting unity for each of the CL ' s and by replacing 
the in's in the expressions for the e's by the magnitude of cjo. This gives: 


4 9* 4 ? • "4 2 * 

D • i (T) <- [cd + 2cn uo + cjo + 4cd“cd + 2co cjo + to + 4crTcD + ao ao ] 

Li 24 


D . (T) < - [3co^ + IOcjo^oo + 4co ao] 

1 24 


( 22 ) 


Similarly (T) and D^ 3 (T) have an upper bound identical to that of equation (22). 
This upper bound gives a means of assessing the update formulas in connection 
with a particular time interval T and mission profile, i.e., relationship 
between no, to, cjo and time. The final error of the C matrix can in principle 
be evaluated by the integral: 
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0 


t f inal 

f [D(t,T)]dt 


where ^- s t ' le e l a P se d mission time. The uncertainty of the spacecraft 

attitude at time t . . , is in turn a function of the final errors in C. 

final 

Of much greater interest than the analytical error expressions derived 
above, are the actually observed errors resulting from the computer simulations. 

3.4 Timing Considerations 

The rationale behind the utilization of a general purpose whole number 
computer to perform the navigation functions in a strapdown navigation system 
is that such a computer must necessarily be included in the spacecraft 
to perform other necessary on-board functions. The argument is made that 
if the percentage of computer time required to perform the navigation function 
is sufficiently small such that the other functions are not adversely affected, 
then the special purpose digital differential analyzer (DDA) which is normally 
associated with the strapdown system can be eliminated from the spacecraft. 
Assuming that a whole number algorithm of updating the cosine matrix is 
sufficiently accurate, the problem reduces to one of comparing an estimate 
of the amount of computer time required to perform the algorithm with the 
amount of excess time capacity of the guidance computer. 

If the Block II AGC as it is presently conceived were required to 
perform the full third order update calculations at the rate of, say, 

10 complete updates per second, then rough estimates indicate that the 
AGC would be saturated performing this task alone. However, it is estimated 
that the AGC could perform an economized version of the third order update 
formula, N^CT), in less than eight milliseconds. (For a complete description 
of the economized form, see Section 5.8.1 Comput er Word Lengt h, page .) 

It is estimated that such an economized form would require less than 
ten percent of the AGC's computing time. A rough estimate of the types 
and number of instructions required by the AGC to perform the third order 
update calculations is given in Table 4. 
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Number 


Function Performed 


Total Memory Cycle Times 


3 

3 

3 

9 

3 

3 

3 

30 

9 

27 

27 

18 


Read 9 . 

l 


9 

Shift Once 


12 

Square 


18 

Cross Multiplies 


54 

Sum of Squares 


6 

9^9 . Terms 

l 


18 

Multiply by 2 


6 

Double Precision 

Adds 

90 

Adds 


27 

Multiplies 


108 

Double Precision 

Adds 

108 

Exchanges 


36 


APPROXIMATE TOTAL 


600 MCT « 7 msec. 


TABLE 4. ROUGH ESTIMATE OF INSTRUCTIONS AND TIMES REQUIRED TO PERFORM 
THIRD-ORDER UPDATE CALCULATIONS IN AGC 
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IV. DIGITAL DIFFERENTIAL ANALYZER (DDA) 


The basic approach to maintaining an updated cosine matrix using 
DDA techniques is to solve equation (4c). Using rectangular rule integration, 
the DDA updates the transformation matrix by solving the difference equations: 


AC. . = C. 
ij i»J+l 


B,j+2 " C i,j+2 B,j + 1 


where 43 . is that angle change of the spacecraft about the j ^ axis 

b , J 

which results in one pulse of a pulse torqued gyro. 


Solution of equation (23) requires that the update cycle time of the DDA 
be sufficiently short such that not more than one 40 change is observed 
by an axis gyro during the update cycle; i.e., 40 pulses about any given 
axis must be processed one at a time and in the order observed. Examination 
of equation (23) indicates that, since addition of angle changes about three 
different axes is non- commutative the updated transformation matrix is 
dependent upon the order in which the individual elements of the matrix 
are updated. This order dependency of the updating procedure of the DDA 
introduces into the updated matrix an inherent inaccuracy which is a function 
of the updating procedure and of the particular flight profile. 

An analysis of various updating procedures for a DDA and of the errors 

(3) 

associated with each of these procedures has been conducted by R. M. Hession , 
the principal results and recommendations resulting from Hession' s analysis 
were utilized in this study as a basis for comparing the performance of a 
whole number updating algorithm with the performance of a DDA. Hession 
concludes that, considering the tradeoffs involved between required accuracy 
and machine speed, the optimum configuration of a DDA is one designated 
as Serial-Parallel (± separately; with reversal rule). Under this organization, 
a complete update of the transformation matrix consists of three partial 
updates; the DDA must thus operate at a cycle time sufficient for the 
three partial updates to be completed between successive 40 changes about 


(3) 


R. M. Hession, R-481, Analy s is o f a Transfo rma tion Computer U sed 
w ith a Gimballess I MU, MIT Instrumentation Laboratory, January, 1965. 
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any axis. To minimize the error resulting from the order in which the 
updates are performed, the update order is reversed after an angular 
change of is observed about any body axis. 

The difference equations to be updated by a Serial-Parallel organized 
DDA are shown below. To simplify the notation, tu is used instead of 
C _ (K + n) refers to the value of the element after having been updated 

n times. The difference equations are: 


c il 

(K 

+ 

1) = 

= C il 

(K) 






C i2 

(K 

+ 

1) = 

= C i2 

(K) ■ 

' °il 

(K) 

h 3 



C i3 

(K 

+ 

1) = 

= C i3 

(K) + C il 

(K) 

h 2 



C il 

(K 

+ 

2) = 

= C il 

(K + 

1) + 

C i2 

(K + 

D h 3 


C i2 

(K 

+ 

2) = 

" C i2 

(K + 

1) 




(24a) 

C i3 

(K 

+ 

2) = 

‘ C i3 

(K + 

1) - 

C i2 

(K + 

1) h 1 



il 

(K + 

3) = C n 

(K + 

2) 

- C i3 

i2 

(K + 

3) - 

(K + 

2) 

+ C i3 

i3 

(K + 

3) - C u 

(K + 

2) 



(K + 2) h 2 
(K + 2) h 1 
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Upon reversal, the equations become: 


C n (K + 4) 
C. 2 (K + 4) 
C- 3 (K + 4) 


c il 

(K + 

3) - C. 3 

(K + 

3) h 

C i2 

(K + 

3) + C. 3 

(K + 

3) h 

C i3 

(K + 

3) 




C. x (K + 5) 
C i2 (K + 5) 

C i3 (K + 5) 


C. . (K + 4) + C.„ (K + 4) h, 
ll J 

C. 2 (K + 4) 

C. 3 (K + 4) - C. 2 (K + 4) h 1 


C 

C 

C 


il 

(K + 

6 > = C il 

(K + 

5) 



i2 

(K + 

6) = C. 2 

(K + 

5 > - c n 

(K + 

5) h 

i3 

(K + 

6) = C. 3 

(K + 

5) + C n 

(K + 

5) h. 


(24b) 


The set of equations (24) were used to describe a DDA subject to the 
following mechanization rules. Each of the elements of the transformation 
matrix consist of two finite length computer words, Y and R. Only the 
Y words were used in the multiplications with the products added into the 
appropriate R register. The lowest order "slot" of the Y word equals the 
magnitude of 1A©. (The terminology "slot" is introduced because the value 
of A© which was utilized, 1/4 milliradian, is not representable by a negative 
integral power of either 10 or of 2. Thus while in most DDA's 149 corresponds 
to the lowest order bit of a binary register, the value of A3 chosen 
for the simulations prohibits the normal use of the word "bit" for the purposes 
of this study.) The R register is restricted in magnitude to be less than 
1 A© ; when the R register exceeds 1A© , an overflow of 149 is affected into 
the corresponding Y register. Under the above form of mechanization, a 
typical update equation from the set of equations (24) becomes: 
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R. 2 (K + 1) = R. 2 (K) - Y ±1 (K) h 3 

If R i2 (K + 1) exceeds 140, then Y^ (K + 1) is incremented by 140 and 
R_ 2 (K + 1) is decremented by 140. 

In addition to the mechanization described above, a roundoff rule 
was employed. Before using a Y word in a multiplication, the corresponding 
R word was checked to determine the value in the R word. If R equalled 
or exceeded 1/2 40 , then Y was incremented by 140 before being used in 
the multiplication; otherwise the value of Y was not altered. In neither 
case was the value of Y modified as it appeared in the Y register. 

All of the DDA results obtained during the simulations resulted from 
the DDA as mechanized above where 140 = 1/4 milliradian. At any given 
instant in time, the value of an element of the transformation matrix is 
equal to the algebraic sum of the contents of the corresponding Y and R 
registers . 


V. COMPUTER SIMULATION 


5.1 Goals of Simulation 

Extensive simulation has been performed on a Honeywell 1800 computer 
to evaluate the algorithm developed in Section 3.1 and especially to 
obtain a more precise feel for the behavior of the accumulated error of 
the C matrix. The simulations were performed in floating point arithmetic 
with a mantissa of 10 decimal digits and an exponent of 2 digits. 

The simulation programs were designed essentially to: 

a. Simulate rotational velocities and accelerations incurred by a 
spacecraft for any specified flight profile. 

b. Determine changes of spacecraft attitude angles (9's) about 

each of the spacecraft's axis for consecutive sampling time intervals 
of length T for the duration of the flight profile. 

c. Update the third order N matrix of Table 3 at time intervals T. 

d. Determine the true C matrix as a function of time by utilizing 
knowledge of the flight profile to solve equation (5). 

e. Determine the error matrix E(t) by comparing the matrices resulting 
from steps c and d. 

It must be emphasized that the simulations were concerned only with determining 
the efficacy of the algorithms as an analytical method of maintaining the 
C matrix. No effort was made to solve the navigation and guidance equations 
of the simulated flight profiles. Thus this study at best represents an 
effort to investigate only one of the many problems associated with the 
strapdown configuration. 

The final step of the study was to simulate the performance of the 
DDA as represented by the set of equations (24) for certain of the profiles 
and to determine the error matrix resulting from the DDA updating technique. 
These simulations permit a comparison of the performance of the DDA with 
the performance of the whole number algorithms. While results of the DDA 
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simulations are included, the remainder of this report is concerned principally 
with a discussion of the whole number algorithm. 

5.2 Determination of True C Matrix 

The principal uncertainty in the results of the many simulated flight 
profiles is the accuracy of the standard solution against which the results 
are compared. Theoretically, the solution to equation (5) would provide 
the precise standard which is desired; in practice, however, the approximations 
introduced by the computer differential equation subroutine make the accuracy 
of the solution questionable. 

The differential equation subroutine utilized in this study permits 
the solution of any set of simultaneous differential equations of the form 
of equation (5) provided that the highest derivative is piecewise continuous 
and that the locations of the discontinuities are known in advance. To 
relate these restrictions to the problem at hand, it is noted that the LEM, 
and in fact the great majority of present day maneuverable spacecraft, is 
attitude controlled by the thrusting of reaction jets. Throughout the 
simulations, the attitude jets were assumed to be capable of existing 
in only one of two states, "on" or "off." When turned "on" the jets provide 
a thrust which results in a constant angular acceleration; when "off" the 
jets provide no angular acceleration. Because the transition between "on" 
and "off" is assumed to occur instantaneously, the angular accelerations 
measured by the spacecraft - and hence indirectly the elements of the matrix 
[C(t)] as updated by the computer - are discontinuous at the time the attitude 
jets are switched. In order to correctly evaluate equation (5) using the 
differential equation subroutine the times of such discontinuities must 
be known in advance. 

It has been assumed in this study that the differential equation 
solution to equation (5) for each flight profile yields an accurate C matrix 
against which the results of the algorithm can be compared. The accuracy 
of the differential equation subroutine utilized depends upon the size of 
an incremental interval of time At, during which the dependent variable 
must be continuous. To determine the validity of the differential equation 
solution for a given flight profile, the subroutine was run several times 
with decreasing values of the increment At; the resulting solutions were 
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checked for convergence. As an example, at the termination of a particular 
40 second mission, the convergence of the solutions corresponding to decreasing 
time increments At is depicted in Table 5. 

The convergence indicated by Table 5 is typical of convergence observed 
for other flight profiles and suggests that a time increment of At = .0015625 second 
provides a sufficiently accurate solution to equation (5). However, it 
was discovered during the simulations that for flight profiles exceeding 
a duration of 100 seconds, the requirements imposed by the small value 
of At exceeded the single precision capabilities of the computer; it was 
also observed that solutions with a time increment of .0015625 second 
required an unreasonable amount of computer time in proportion to the scope 
of this study. Hence for all flight profiles the "true" C matrix was obtained 
by solving the differential equation (equation (5)) with an incremental 
time interval At = .003125. The resulting C matrix can be considered accurate 
to at least the sixth decimal place. 

5.3 Characteristics of Simulated Flight Profiles 

The flight profiles which were simulated in this study fall into two 
basic categories: 

1. Missions which in the most general case consist of alternate polarity 
acceleration pulses applied independently to the attitude jets of 
each of the three spacecraft axes. 

2. Profiles which represent a typical LEM mission. 

The following constraints were imposed upon the spacecraft maneuvers called 
for in the simulations: 

1. All angular accelerations about each axis were of constant magnitude, 

3/4 radian per second per second. 

2. For profiles of the first category, rotational velocities about 
each axis were limited to magnitudes of 20° per second or less; 
maneuvers in the LEM missions were limited to rotational velocities 
whose magnitudes were 10° per second or less. 
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TABLE 5. CONVERGENCE AT 40 SECONDS OF DIFFERENTIAL EQUATION SOLUTION OF 
COSINE MATRIX FOR DECREASING TIME INCREMENTS, At (At MEASURED IN SECONDS) 
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Characteristics of successive runs were dictated by the desire to isolate 
the effects of the various parameters which influence the capabilities of 
the basic algorithm. 

5.4 Effect of Vehicle Rotation 

The characteristics of the first few simulations were designed to 
isolate the effects of vehicle rotation upon the matrix [E (T) ] given in 
Table 3. A single acceleration pulse, commencing at time t = 0 and of 
a specified duration, was applied to the x axis attitude control jets; 
subsequent to the termination of acceleration, the spacecraft was allowed 
to rotate freely at a constant angular velocity for a duration of 200 seconds 
during which time the third order N matrix of Table 3 was updated at sampling 
time intervals of 0.1 second. Periodically during the 200 seconds, the 
updated N matrix was compared with the true transformation matrix [C ( t ) ] 
(differential equation solution) and the resulting error matrix 
[E(t)J = [C(t)] - [N(t)J was determined. The magnitude of the elements 
of [E(t)] represents the degree to which the updated N matrix approximates 
the M matrix of equation (6). 

The profile described above was simulated for the acceleration pulse 
lasting .025 second, 0.25 second and .465625 second. (The last value 
represents the approximate time which it would take a body under an angular 
acceleration of 3/4 radian per second per second to achieve a rotational 
velocity of 20° per second.) The behavior of one element, e 22’ t ^ le resu lti n g 
error matrix for each of the three profiles is shown in Figure 1. It 
should be noted that there is nothing unique about the element e 22 ! i- ts 
behavior is simply typical of - but not identical to - the other elements 
of the error matrix. 

There are several interesting properties about the functions shown 
in Figure 1. According to the error matrix of Table 3, it was expected 
that the element should be essentially a constant for the duration 

T 4 

1 4 

of the profile since a 99 of Table 3 reduces in this case to - co . It 

24 X 

is apparent from Figure 1 that the truncated error matrix of Table 3 does 
not adequately represent the behavior of the basic algorithm; evidently 
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higher order terms must be considered. The periods of oscillation of the 
error functions are identical to the time it takes for the spacecraft to 
rotate 360° while the peaks of the error function oscillations grow in a 
linear fashion which is apparently a function of the speed of rotation. 

The unfortunate result is that the error peaks appear to be monotonically 
increasing. From the three functions depicted, no simple relationship 
between the rotational velocity and the growth of the error peaks has been 
determined. One is tempted to conclude from Figure 1 that the errors 
resulting from the basic algorithm are functions of the spacecraft velocity 
and attitude. It should be remembered, however, that the differential 
equation solution was shown to converge only to the sixth decimal place 
for At = .003125; therefore, considering the magnitude of the error, one 
might question the accuracy of these initial conclusions. 

To ensure that simultaneous rotation about each of the body axes does 
not adversely affect the performance of the basic algorithm, a flight profile 
similar to the above was simulated. This profile consisted of applying 
a .025 second acceleration pulse about each of the body axes at time t = 0 
and then allowing the spacecraft to rotate freely for 200 seconds. The 
behavior of element the resulting error matrix is shown in Figure 2. 

While unquestionable conclusions cannot be drawn from one simulation, 
nevertheless comparison of Figure 2 with Figure la indicates that simultaneous 
rotation about the three body axes does not significantly affect the accuracy 
of the basic algorithm. 

5.5 The Basic Profile 

From the results observed for constant rotation of the spacecraft, 
it became apparent that more sophisticated maneuvers must be studied. 

A flight profile was designed which consisted essentially of limit cycle 
maneuvers performed about each of the body axes; because this profile was 
the basis of the great majority of the simulations, it will hereafter be 
referred to as the basic profile. The characteristics of the accelerations 
applied about each of the body axis are shown in Figure 3. The sequence 
of the pulses about the z axis permits acceleration to the maximum allowable 
rotational speed of 20° per second, free rotation at this velocity for 
a period of time, followed by deceleration to zero rotation about the 
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z axis. The vehicle is assumed to be rotating at time t = 0 about the x and 

the v axes with a velocity oo = CD = -.009375 radian per second; thus the 

x y 

x and y axis limit cycle maneuvers are centered about the respective axis. 

The basic algorithm was used to update the N matrix during simulations 
of the above described profile for elapsed mission times of 1000 seconds. 
Sampling time intervals of 0.1, 0.05, and 0.025 second were employed. A 
crude graph of the behavior of one element, e^, of the error matrix 
[E(t)] = tC(t)] - [N(t)] for each of the three values of the sampling 
time interval is shown in Figure 4. Figure 5 is a more detailed presentation 
of the behavior of e^ for the first 200 seconds of the mission with a 
sampling time interval T = 0.1 second. 

Figure 4 substantiates the prediction that a reduction in size of the 
sampling time interval results in a corresponding decrease in the magnitude 
of the errors. According to the error matrix of Table 3, a reduction in 
size of the sampling time interval by a factor of two should reduce the 
error by a factor of sixteen. While it is not immediately apparent from 
the functions of Figure 4, comparison of corresponding raw data points 
indicates that halving the sampling time interval results in a reduction 
of error magnitude by a factor of only five. This result tends to substantiate 
the earlier conclusion that at least some of the higher order effects which 
were omitted in the basic algorithm and the error matrix are not truly 
negligible . 

There are several interesting observations which can be drawn from the 
function shown in Figure 5, which is an expansion of the first 200 seconds 
of the profile. The frequency of the sinusoidal type pulses, which evidently 
result from high speed rotation about the z axis, is identical to that 
observed in Figure lc. Since the magnitude of rotatation is the same 
for both cases, this is not an unexpected result; however the magnitude 
of the error now indicates that the differential equation routine is not 
the cause of oscillatory behavior of the error function. Where the peaks 
of the error in Figure lc grow linearly, such is not the case for the basic 
profile. In fact, according to Figure 4, the magnitude of the error peaks 
is well bounded. The most significant result which is apparent in Figure 5 
is the difference in the order of magnitude between the errors observed 
for this profile and the error function of Figure lc even though the magnitude 
of body rotation is the same for both profiles. Such a discrepancy can only 
be caused by one or more of the following: 
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1. Simultaneous acceleration about one body axis and rotation about 

the other body axes; i.e., occurrrence of non-synchronous accelerations 
about the three body axes. 

2. Occurrence of limit cycling or repetitive accelerations. 

It should be remembered that simultaneous rotation of small magnitude 

about the three body axes did not previously appear to affect the performance 

of the basic algorithm. 

To determine the effects of the limit cycle maneuver, a profile similar 
to the basic profile was s imululated ; in this profile, however, no accelerations 
were applied to the z axis. The error function e^^(t) which results from 
limit cycle maneuvers about the x and z axes is shown in Figure 6. Although 
the speed of rotation about the two axes of the spacecraft was identical 
to that of the earlier profile, the magnitude of the errors is nevertheless 
considerably larger than that observed in Figure la. Because the x and y axis 
accelerations are synchronous, the error shown in Figure 6 can only be due 
to the repetition of accelerations. It appears that the errors are somewhat 
cumulative; however the periodicity of groups of three error pulses remains 
unexplained at this time. 

Further evidence that the performance of the basic algorithm is 
influenced by the occurrence of repetitive accelerations is presented in 
Figure 7 which shows the error function e^(t) resulting from simulation 
of the z axis accelerations only of the basic profile. For the profile 
in which one .465625 second acceleration pulse was applied to the z axis 
attitude jets followed by a 200 second period of free rotation (profile 
for Figure lc) , the error function e^Ct) was identically zero. Thus 
the existence of the error function shown in Figure 7 can be caused only 
by the repetitive accelerations. 

5.6 Non-Synchronous Accelerations 

The simulations considered thus far have consisted of applying accelerations 
simultaneously to various combinations of the spacecraft axes. To investigate 
the effect of asynchronous accelerations, two profiles were simulated. 

The profiles consisted of accelerations of alternating polarity applied 
to the x, y, and z axes at multiples of 4, 5, and 7 seconds respectively; 
in one profile the accelerations lasted for .025 second while in the other 
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profile the pulses were of .465625 second duration. The error functions e 31 (t) 
for the .025 and the .465625 second profiles are shown in Figures 8a and 
8b respectively. The functions shown in Figure 8 substantiate the previous 
conclusion that the magnitude of error is dependent in part upon the speed 
of rotation of the spacecraft. It is significant to note, however, that 
the error shown in Figure 8b is of smaller magnitude than the corresponding 
error of the basic profile, even though the profile for Figure 8b calls 
for high speed rotation about each of the spacecraft axes while the basic 
profile calls for high rotation about only the z axis. Such a result was 
completely unexpected and at the time unexplained. 

5.7 Extensions of the Basic Algorithm 

Of particular concern to the design of the LEM navigation and guidance 
system is the magnitude of errors resulting from seemingly non-violent 

maneuvers. Since the updated N matrix is an approximation to a matrix 

_2 

of direction cosines, an error of 3 x 10 (see Figure 4) can represent a 
spacecraft attitude error of almost 2 degrees, an error which is unacceptable 
for the LEM mission. In at attempt to reduce the size of the errors while 
simultaneously gaining more insight into the characteristics of the basic 
algorithm, several extensions of the algorithm were developed and simulated. 
These modifications and the results of their simulation are described below. 

5.7.1 Reduction of Sampling Time Interval 

The error matrix of Table 3 predicts, and the error functions of 
Figure 4 verify, that a reduction of the sampling time interval results 
in a reduction of the error magnitude. However, if the sampling time interval 
is made small enough to ensure that the updated N matrix closely approximates 
the true solution, an unreasonable computational load is placed upon the 
navigation computer. An attempt was made to realize the reduction of error 
magnitude by sampling attitude angle changes at relatively short time intervals 
while performing update calculations only at longer time intervals. The 
update calculations are of course more complex than the elements of the 
N matrix shown in Table 3. 
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Referring to Figure 9, let the long update time interval T be divided 
into two shorter time intervals, each of length T/2. 



Fig. 9. Reduction of Sampling Time Interval 

Within each interval T, denote the change of attitude angle about the i t ' 1 body 

axis during the first interval T/2 by Oh and the change of angle during 

the second interval T/2 by 3^; thus the total angle change, 9^, equals 

OL. +3.. If the N matrix of Table 3 is updated at intervals T/2 rather 
r a. 

than intervals T, the N matrix at time T is: 

[N(a, 3, T)] = [N(a, T/2)] [N(3, T/2)] (25) 

where [N(a, T/2)] and [N(6, T/2)] have the same elements as given in Table 3 
for [N(9, T) ] except that 9 becomes a and 3 respectively. By sampling 
the angles a and 3 at times T/2 the result given in equation (25) can be 
obtained by updating a new N matrix, [N 1 ], periodically at time intervals T. 

The expression for [N 1 ] is of course: 

[N'(a, 3, T)] = [N(a, T/2)] [N(3, T/2)] (26) 

and includes terms of the sixth order rather than the third order as given 
in Table 3. (The expansion [N(QI, T/2 )][n( 3, T/2)] is rather tedious and, 
since it is of no real significance, is not included in this report.) Thus 
the effect of updating the N matrix at shorter time intervals can be realized 
by measuring the angle changes at shorter time intervals T/2 while performing 
update calculations periodically at longer time intervals T. 

The above procedure was simulated with the following exception: in 
order to keep the computational load on the navigation computer at a reasonable 
level, it was decided not to include in the elements of the N' matrix terms 
which were of the third order or greater. The results of the simulation 
of the N' matrix approximation are not plotted but in general the elements 
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of the error matrix resulting from this procedure are somewhat smaller 
than those resulting from the basic N matrix approximation; however, the 
reduction of the error matrix is less than one order of magnitude. Considering 
that the computational load is higher for the N' matrix than for the N matrix, 
such results are not encouraging. 

From the expansion of the product [N(a, T/2 )][n( 8, T/2)], it was observed 
that the leading terms of the difference ([N(a, T/2 )][n( 8, T/2)]) - [N(9, T)] 
(remembering that 9 = a + 8) were terms typically of the form (<X 8. - Ct. 8.)/2. 
Since these terms do not represent an unreasonable amount of computation, 
a logical suggestion is to include these terms in the elements of the third 
order N matrix of Table 3 to determine whether or not they contribute significantly 
to the reduction of error observed in Figure 4 when angle changes are sampled 
at twice the update rate. The results obtained from the simulation of this 
amended N matrix approximation indicate that the elements of the resulting 
error matrix are slightly smaller than the error terms resulting from the 
N' matrix approximation. Again, however, the decrease in error magnitude 
is less than one order of magnitude. 


5.7.2 Fourth Order N Matrix 

The error matrix shown in Table 3 results from truncating the expansion 
of the N matrix elements of the third order. To verify that the updated 
third order N matrix is in fact a reasonable approximation to the cosine 
matrix, an N matrix was constructed whose elements include the fourth order 
terms necessary to eliminate the fourth order terms of the error matrix. 

t 

The fourth order N matrix was then used in simulated profiles to approximate 
the cosine matrix. That it is unnecessary to include fourth order terms 
in the N matrix was demonstrated by the fact that the resulting errors 
were at least 907o as great as the errors observed for the third order N matrix. 

The results obtained for the fourth order N matrix in turn suggest 
that the necessity of including third order terms in the N matrix is questionable. 
Profiles were simulated, wherein various combinations of the third order 
terms of the N matrix of Table 3 were not included in the elements of 
the update matrix; in one simulation, no third order terms were included. 

The results of these simulations showed that the omission of the third 
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order terms of the basic algorithm results in error functions which are 
one order of magnitude greater than the error functions resulting from a 
full third order N matrix. It thus appears that, at least for the flight 
profiles which were simulated, an updated full third order N matrix represents 
a reasonable approximation to the direction cosine matrix with a corresponding 
acceptable level of computational complexity. 

5.7.3 Interrupted Sampling Time Interval 

The derivation of the basic algorithm results from the Taylor series 
expansion of the changes of spacecraft attitude angles during a specified 
time interval as given by equation (13) . For equation (13) to be a valid 
representation of 9, the function 9(t) and all of its time derivatives 
must be continuous throughout the time interval T. However, for the method 
described for implementing the basic algorithm, the requirement for continuous 
derivatives is not necessarily met. Permitting the spacecraft attitude 
control jets - and hence the body angular acceleration - to be in only 
one of two states, on or off, introduces a discontinuity in the angular 
accelerations at the time the control jets are switched. Unless these 
discontinuities occur at the initiation of a sampling time interval, the 
function 9(t) has discontinuous derivatives ard the expansion of equation (13) 
over the sampling interval becomes invalid. 

A necessary condition to ensure that the function 9(t) has no discontinuous 
derivatives during a sampling interval is that changes in spacecraft acceleration 
be made coincident in time with the beginning of a sampling time interval. 

For any realistic flight profile, it is impractical to predict the exact 
times when changes in acceleration will occur; it is therefore impossible 
to determine a priori a fixed value T for the sampling time interval which 
guarantees coincidence between the sampling interval and acceleration change 
for the entire mission. The obvious solution is to force, at the time of 
acceleration change, the current sampling interval to be terminated and 
the subsequent interval initiated. 

Several practical methods of implementing the interruption of the 
sampling interval can be suggested. If the navigation computer is also 
performing the guidance functions of the spacecraft, knowledge of the 
attitude jet firing is implied; otherwise an interrupt signal from the 
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guidance equipment to the navigation computer is required. In either case, 
the navigation computer can terminate the current sampling time interval, 
perform update calculations and initiate the subsequent sampling time 
interval; thus no discontinuous derivatives are permitted to occur during 
a sampling time interval. On the other hand, in any realistic space mission, 
there are bound to be rotational accelerations of the spacecraft which 
will not be initiated or observed by the on-board computer. Fuel slosh, 
motion by the spacecraft occupants, etc., result in changes in rotational 
accelerations about which the guidance computer has no knowledge. Thus 
complete coincidence between sampling time intervals and changes in vehicle 
accelerations cannot be simply assured. It might be noted, however, that 
these sources of accelerations are not as sharp as the jets and presumably 
will not introduce large additional errors. 

Although only one method of interruption was simulated, it is felt 

on the basis of earlier simulations that other methods would yield essentially 

the same results. Under the simulated procedure, the length of the sampling 

time interval is set at some constant value T and, in the absence of 

acceleration changes, update calculations are performed periodically at 

time intervals T as usual. When an acceleration change occurs at time t 

during the n^ sampling time interval (n = integer) ( the n^ interval is 

s t 

terminated at time t and the length of the n + 1 interval is set at nT - t. 

Until the next acceleration change occurs, all subsequent intervals are 

s t 

of length T. (The length of the n + 1 interval could have been set at T.) 
Changes in spacecraft attitude angles are noted and update calculations 
are performed periodically every nT seconds for the duration of the flight 
profile and in addition are performed at every change of acceleration. 

Although the interruption of the sampling time interval assures that 
the M matrix of equation (6) can in theory be approximated by an N matrix 
such as that given in Table 3, nevertheless the interruption introduces 
an additional error which results from the fact that all sampling intervals 
are not of the same length. To illustrate this inherent error, consider 
the case where the n t ^ sampling interval is of length T while the n + 1 st 
interval is of length t. For the sake of notational convenience, let time 0 
represent the termination of the n^ interval and the initiation of the n + 1 st 
interval; in other words, at time t a change in acceleration occurs resulting 
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in two consecutive sampling intervals of different lengths. The time 

sequence is illustrated in Figure 10. It is convenient now to investigate 

some of the functions of Table 2 as they are calculated at the end of the 
s t 

n + 1 interval. For example: 
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2 2 
Tt tT 

= - + - [to. (0)o>. (0) - to. (0)oo. (0) ] (27d) 

2 2 1 J 

Inspection shows that equation (27d) is identical to the corresponding 
function given in Table 2 only when t = T, a condition which is not possible 
when an interrupt occurs. Thus the interrupt inherently introduces an error 
not shown in Table 3 which in this case is represented by the factor 

3 2 2 

T - - (Tt + tT ) . In general the error is a function of the difference 

2 

in length between consecutive sampling time intervals. 

Updating the full third order N matrix given in Table 3 with an 
interrupted sampling time interval was simulated using the basic profile. 

The normal sampling time interval was 0.1 second, the same value used in 
the earlier simulations. The error functions e^Ct) resulting from the 
interrupted simulation is shown in Figure 11. The principle result to 
be noticed from a comparison of Figure 11 with Figure 5 is that the error 
has been reduced by three orders of magnitude. Such significant reduction 
is not limited to this profile; comparison of the interrupted update method 
with the basic algorithm for other profiles revealed similar results. 
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5.8 


Practical Considerations 


5.8.1 Computer Word Length 


A full third order update N matrix approximation evaluated ten times 
a second and at every change of acceleration results in an error matrix 
which corresponds to at least milliradian accuracy for the profiles which 
were simulated and in considerably better accuracy for the LEM profile 
which was simulated. To achieve such accuracy with the 15 bit word length 
of the Apollo Guidance Computer, calculations would have to be performed 
in double precision form. A rough estimate indicates that if the interpretive 
mode were employed, the computer would be fully occupied with this job alone. 

In order to reduce the time requirement levied against the AGC by the update 
formula, the calculations are expressed in a more convenient form. 

The variables are first scaled such that the number of double precision 

additions is minimized. Terms of the N matrix in 9 occupy the higher component 

2 2 

of the double precison word; terms in 9^, 9^9^. and ® occupy the lower 
component. The maximum value of the 0^ is scaled by choosing the interval T 
in accordance with the maximum angular velocity of the spacecraft. These 
terms are accumulated to form the nine elements of the matrix [N-l] where 
[i] is the identity matrix. The matrix multiplication [C] x [N-l] is performed 
in single precision using the upper components of each matrix and the resulting 
double length product is added to [c]. The lower component of [N-l] is saved 
to be added to the next sampling period's [N-l]. The process is like that 
of the DDA where the lower component of a double precision word is saved 
and accumulated at each step. Because this form of mechanizing the computational 
procedure yields results which are less accurate than those resulting from 
precision, it will hereafter be referred to as computation of 1 1/2 precision. 

The 1 1/2 precision form of computation described above was implemented 
on the Honeywell 1800 with two slight variations. 


1. The 15 bit word length of the AGC was simulated as a five digit 

decimal word length. Thus the accuracy of the simulated solutions 
should be less than the true solutions for the described procedure. 
(Because the N matrix is a matrix of cosines, the first digit of 
either the binary or decimal word must have a magnitude of either 
0 or 1 leaving 4 decimal digits and 14 bits to provide the precision 
of the cosine.) 
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2. All terms of the N matrix were accumulated in the lower component 

of the double precision word but were allowed to overflow into 

the higher component of the double precision word. However, 

because of the floating point arithmetic employed by the simulation 

computer, the higher component of the simulated double precision 

word always contains five significant digits of the cosine approximation. 

Thus the simulated procedure is accurate to the extent that terms 
2 2 

of 6^ and 9 0^ do not contribute to the first five significant 
digits of the cosine approximation. The effects of the 1 1/2 
computational precision will be discussed with the results of the 
simulated LEM mission. 

5.8.2 Gyro Limitations 

The results presented to this point are predicated upon the assumption 
that the changes in the spacecraft attitude angles during each sampling time 
interval are known precisely, i.e., that the gyros by which the angle changes 
are measured provide a continuous readout of angle change data to the navigation 
computer. In practice, of course, this situation is not realized. In the 
case of the LEM mission, the gyros are pulse torqued gyros which require 
one output pulse from the computer for each change of attitude angle £6, 
a positive pulse for a net positive angle change and a negative pulse for 
a net negative angle change. Because angle changes are algebraically accumulated, 
are measured with respect to a fixed reference, and can only be measured 
to the nearest A© through which the vehicle has rotated, it is possible for 
the spacecraft to rotate between the angles +A0 and -A3 with no pulses being 
applied to the gyros. This range of angles is known as the dead band of 
the gyros and introduces a memory type effect into the determination of 
angle changes. A dead band of 1/2 milliradian (AS = ± 1/4 mr.) was introduced 
into the simulations of all profiles; the effect of introducing imperfect 
gyros was observed to be essentially independent of the mission profile 
and is therefore discussed only with the results of the simulated LEM mission. 
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5.9 LEM Profile 


The fundamental purpose of the entire investigation was to determine 
the practicality of utilizing the N matrix given in Table 3 to maintain 
spacecraft attitude during the LEM mission. The LEM mission profile which 
was simulated is in fact a fairly simple approximation to an actual LEM 
profile, the simplicity being a result of the author's ignorance of the 
detailed flight profile. The approximation will, however, suffice for 
the purposes of evaluating the results of the update approximation. 

The simulated LEM mission profile consisted of maneuvering the spacecraft 
about the pitch (x) axis while limit cycling the LEM about the roll and yaw 
(y and z) axes. The maneuvers performed about the pitch axis can best be 
described with reference to Figure 12a and consist of the following: 

1. At time t = 0, the LEM leaves the orbiting platform with an inertial 
pitch rate of 0.1°/second, allowing the LEM to retain local orientation 
with the moon. 

2. At t = 300 seconds, the LEM pitches 20° in preparation for approaching 

the moon; this attitude is held for 150 seconds while the LEM 

descends towards the surface. The pitch of 20° is made at the 

maximum allowable rotation of 10°/second and maximum accelerations 
2 

of 3/4 radian/second . 

3. At 450 seconds, the LEM pitches 60°, at which point it is oriented 
with the local vertical of the moon. 

4. For 120 seconds, the LEM hovers in a vertical attitude over the 
landing spot. 

5. Failing to find an acceptable landing point, the LEM aborts a 
lunar touchdown, pitches 60° and lifts off to a rendezvous with 
the orbiting platform; during this period of the mission, a pitch 
rate of 0.1° is again employed to maintain local orientation. 
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Limit Cycle Maneuvers for Simulated 






6. The mission terminates at t = 1000 seconds at which time rendezvous 
with the orbiting command module occurs. 

7. Throughout the mission, limit cycling as shown in Figure 13b is 
occurring about the roll and yaw axes of the LEM. 

As a result of knowledge gained from the earlier simulations, maintenance 
of the LEM attitude for LEM profile was simulated using the basic algorithm 
both with and without interrupted sampling time intervals. For each of 
these methods of updating the N matrix, the following combinations of computational 
precision and gyroscope performance were simulated: 

1. full precision, ideal gyros 

2. full precision, gyro readout quantized at 1/4 milliradian 

3. 1 1/2 precision, ideal gyros 

4. 1 1/2 precision, gyro readout quantized at 1/4 milliradian 

It should be mentioned that the above combinations were also simulated 
for the profiles described earlier in this report and the results described 
below were observed for all profiles. 

For the four simulations using the basic update formula, the error 
functions corresponding to the 9 elements of the error matrix are shown 
in Figure 13. Only very general and almost insignificant statements can 
be made about these error functions. The first such statement is that 
the errors resulting from computing in 1 1/2 precision are of the same 
magnitude as those resulting from the utilization of quantized gyros. 

One interesting result is that the errors resulting from the combination 
of computing 1 1/2 precision and of using quantized gyros are not significantly 
greater than the errors caused by either of these two factors separately. 
Furthermore, at least for the simulated LEM mission, the basic algorithm 
utilizing full computational precision and employing ideal gyros yields 
error functions which are of the same order of magnitude as the errors 
resulting from the utilization of 1 1/2 precision and quantized gyros. 

Figure 14 shows the 9 error functions resulting from the simulation 
of the LEM mission using the third order N matrix with interrupted sampling 
intervals. It is significant to note that the errors resulting from the 
utilization of this method with full computational precision and perfect 
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gyros are orders of magnitude less than errors resulting from the other 
combinations of computational precision and gyro readout. Thus the limiting 
factor of the accuracy of the cosine matrix approximation for this profile 
is not the update method itself but rather can be attributed to external 
sources. It is again noted that the effect of 1 1/2 precision is essentially 
the same as the effect of quantized gyros and that these effected do not 
appear to be cumulative. 

To more readily compare the results of the basic algorithm with those 
of the interrupted sampling interval update method, certain of the functions 
of Figures 13 and 14 are superimposed and presented in a common co-ordinate 
system in Figure 15. Figure 15 substantiates the major conclusion which 
was evidenced earlier: elimination of discontinuities in the angular 
accelerations (and the higher derivatives) results in significant improvement 
in the performance of the basic algorithm. 

Figures 13, 14 and 15 reveal the performance of the basic algorithm 
for a mission which somewhat approximates one possible LEM mission. In 
order to delineate the relative performance between the whole number algorithm 
as herein implemented on a general purpose computer and the specialized 
techniques employed by what is felt to be a better-than-average configuration 
of a DDA, the performance of the DDA described in Section IV was simulated. 
The nine elements of the error matrix, the difference between the true trans- 
formation matrix and the transformation matrix as updated by the DDA for 
the LEM profile, are shown in Figure 16. Comparison of Figure 16 with 
Figures 13, 14 and 15 indicates that in general the elements of the trans- 
formation matrix as updated by the DDA are somewhat more accurate than those 
of the whole number algorithm as simulated. 

5.10 Comparative Data for the Basic Profile 

The error functions shown in Figures 13, 14, 15 and 16 present a 
reasonably complete picture of the capabilities of the whole number algorithm 
for the simulated LEM mission. To forestall the possibility of doubt arising 
about the general validity of the information depicted because of the mild 
characteristics of the LEM mission profile, a complete set of error functions 
was obtained for the first 200 seconds of the basic profile. Of all the 
mission profiles which were simulated, this profile resulted in the largest 
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magnitudes of the error functions. However, the maneuvers encountered 
in the basic profile are far more extreme and violent than one would anticipate 
for an actual spacecraft mission; therefore the large errors resulting 
from this profile should not be viewed with excessive alarm. 

Figure 17 shows the 9 error functions of the basic profile resulting 
from three mechanization combinations of the noninterrupted update algorithm. 
Figure 18 shows the 9 error functions of the basic profile resulting from 
simulation of the interrupted update procedure. Figure 19 shows the 9 error 
functions of the basic profile which result from simulation of the DDA 
techniques. The combined set of Figures 13 - 19 present a reasonably complete 
picture of the absolute and relative capabilities of the whole number 
algorithm derived in Section III. It is felt that they, in conjunction 
with the results previously described, provide an objective basis from 
which designers of data processors associated with strapdown navigation 
systems are free to draw their own conclusions. 
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VI. MATRIX ORTHOGANALITY 


The updated N matrix, regardless of how derived, is an approximation 
to the true cosine matrix C given in equation (1), a matrix whose elements 
are the direction cosines of the angles between the axes of the body co- 
ordinate system and the axes of the fixed or inertial coordinate system. 

Since the two coordinate systems are each orthogonal, it follows that the 
matrix C represents an orthorgonal transformation and is therefore orthogonal. 
Furthermore if the C matrix is orthogonal, it follows from equation (6) 
that the M matrix must also be orthogonal. The update formula given in 
Table 3 yields an N matrix which is an approximation to the M matrix; therefore 
[N] should also be orthogonal. 

Consider the matrix: 

[Z] = [N][N] T (28) 

T 

where [N] is the transpose of [N]. If [N] is orthogonal, [Z] is the 
identity matrix; the variation of [z] from the identity matrix provides 
a measure of the degree of orthogonality of [n]. 

It was initially suggested that the Z matrix might provide an evaluation 
of the update formula and might further be used to "correct" the N matrix 
on a real-time basis. In practice, however, it was found that the difference 
between [z] and the identity matrix provides at best a crude indication 
of the effectiveness of the update methods. Examination of the Z matrices 
given in Table 6 reveals the difficulty in constructively utilizing the 
property of orthogonality to correct "the update formulas." For example, 
according to Figure 15, the use of interrupted sampling time intervals 
yields consistently more accurate results than those realized from the 
noninterrupted update technique; yet at the termination of the LEM mission, 
the Z matrix for the noninterrupted algorithm is closer to the identiy 
matrix than is the Z matrix for the interrupted updating. Similarly, according 
to the Z matrices d, e, and f of Table 6 the N matrix of the noninterrupted 
update method is, at the time of maximum error magnitude, less orthogonal 
than the N matrix of the interrupted N matrix at the same time. However, 
at the time of peak error magnitude, the N matrix for noninterrupted update 
formulas is more orthogonal than the N matrix of the interrupted formulas 
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at the time of its largest error functions. Such results indicate that the 
property of orthogonality cannot be constructively utilized. 
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(a) 


Z Matrix at End of LEM Mission Resulting from Basic Update Formula, 
Perfect Gyros and Full Computational Precision. 


1.0000119 
- .00000012 
.00000016 


- .00000012 
1.0000114 
.00000018 


.00000016 

.00000018 

1.000011 


(b) 


Z Matrix at End of LEM Mission Resulting from Interrupted Update Formula, 
Perfect Gyros and Full Computational Precision. 


1.0000122 

.00000044 

.00000045 


.00000044 

1.0000133 

.00000016 


.00000045 

.00000016 

1.000013 


(c) 

Z Matrix at End of LEM Mission Resulting from Interrupted Update Formula, 
Gyros Quantized at 1/4 Milliradian, 1 1/2 Precision. 

"".99955 .000070 .000015 

.000070 .99984 .000084 

.000015 .000084 .99981 


(d) 


Z Matrix at t = 190 of 
Perfect Gyros and Full 

.9998908 

.00001186 

.00002645 


Profile #1 Resulting from 
Computational Precision. 

.00001186 
.9998949 
.00002999 


Basic Update Formula, 

. 00002 645~ 

.00002999 

.99997939 


TABLE 6. TYPICAL Z MATRICES FOR VARIOUS PROFILES 
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(e) 


Z Matrix at t = 190 of Profile #1 Resulting from Interrupted Update Formula, 
Perfect Gyros and Full Computational Precision. 


.9998918 

.00001155 

.00002580 


.00001155 

.9998965 

.00002995 


.00002580 

.00002995 

.99997979 


(f) 


Z Matrix at t = 193 of Profile #1 Resulting from Interrupted Update Formula, 
Perfect Gyros and Full Computational Precision. 


.9998889 

.00001248 

.00002684 


.00001248 

.9998943 

.00003142 


.00002684 

.00003142 

.9999778 


TABLE 6 (cont . ) 



CONCLUSIONS 


Extensive computer simulations have verified that the transformation 
matrix required for attitude reference in a strapdown inertial navigation 
system can for certain missions be updated at relatively long time intervals 
with sufficient accuracy by an on-board general purpose whole number computer. 

The accuracy of the update formulas is strongly dependent upon the characteristics 
of the particular flight profile, specifically upon the spacecraft rotational 
velocity, rotational acceleration, the number of times angular accelerations 
are encountered, and the times at which the accelerations occur. The third 
order update formulas offer a reasonable compromise between computational 
complexity and accuracy of the updated matrix; little improvement is realized 
by using fourth order formulas while significant degradation results from 
second order expressions. Sampling time intervals of the order of 0.1 second 
are sufficiently small to yield maaningfull results; smaller intervals will 
of course yield more accurate results but will also place an increasing 
computational load upon the navigation computer. Real time knowledge of 
the occurrence of discontinuities in the time derivatives of the angular 
velocities can be used to significantly improve the performance of the basic 
algorithm. A sufficient increase in computer word length such that computations 
can be performed in single precision without excessive loss in computational 
accuracy results in improvement in the accuracy of the whole number algorithm. 
Similarly, an increase in precision of readout from the strapped down gyros 
results in improved accuracy of the transformation matrix. 

Relating the results of the study to the AGC Block II computer, it 
appears that the Block II computer yields results in terms of accuracy which 
are barely acceptable. Second generation airborne computers will, however, 
operate at speeds five to ten times faster than the Block II AGC. The increased 
accuracy resulting from shorter sampling intervals thus makes a strapdown 
navigation system employing a general purpose computer very attractive. 
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